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Abstract.  The  ability  to  accurately  simulate  rare  high  energy  collisions  such  as  those  found  in  ionization  and  combustion 
reactions  is  important  to  the  advancement  of  the  study  of  non-equilibrium  behavior  of  these  phenomena.  A  major  difficulty 
for  modeling  these  systems  with  particle  codes  is  that  chain-branching  reactions  result  in  exponential  growth  regimes  for 
certain  species.  This  work  considers  the  use  of  fractional  collision  models  in  place  of  statistically  sampled  collision  models 
combined  with  merging  of  computational  particles  to  reduce  statistical  noise  while  avoiding  runaway  computational  cost. 
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INTRODUCTION 

The  advent  of  the  ’no-time-counter’  (NTC)  collisional  sampling  technique  by  Bird[l]  greatly  contributed  to  the 
proliferation  of  DSMC  codes  used  throughout  the  rarefied  gas  dynamics  field  by  making  physically  relevant  gas 
dynamics  problems  computationally  accessible.  This  was  achieved  by  sampling  only  a  fraction  of  possible  collision 
pairs  and  adjusting  the  probability  of  collision  to  account  for  the  skipped  collisions.  The  technique  produces  physically 
accurate  collision  rates  for  gases  and  gas  mixtures  in  both  equilibrium  and  non-equilibrium [2]  flows.  However  for  some 
types  of  rare  collisions  such  as  inelastic  ionization  processes,  this  approach  results  in  significant  noise. 

For  any  particle  technique,  resolving  the  tails  of  the  velocity  distribution  function  (VDF)  is  necessarily  difficult. 
By  definition,  the  tail  represents  a  large  volume  of  velocity  space  yet  only  a  small  fraction  of  the  total  number  of 
particles.  Many  of  the  special  inelastic  collisions  found  in  plasma  and  combustion  processes  require  overcoming  an 
energy  barrier.  They  therefore  disproportionately  rely  on  collisions  by  particles  in  this  tail. 

When  only  a  few  computational  particles  can  be  devoted  to  resolving  the  tails,  what  should  be  represented  as  a  large 
number  of  low  probability  events  is  instead  represented  by  a  single  infrequent  event  involving  a  large  fraction  of  the 
probability  in  that  region  of  the  VDF.  Suppressing  this  noise  requires  large  numbers  of  particles  so  that  more  than  just 
a  few  particles  resolve  the  relevant  tail  of  the  VDF. 

An  example  of  this  behavior  can  be  seen  in  Bird’s  OD  study  of  ignition  delay  time  using  the  quantum-kinetic 
chemistry  model [3].  The  delay  time  appears  to  converge  only  when  using  greater  than  between  10e5-10e6  particles. 
This  is  a  result  of  the  relatively  low  probability  of  chain  initiation  reactions.  Such  a  high  particle  count  quickly  becomes 
intractable  for  any  higher  dimensionality  flow. 


FRACTIONAL  COLLISIONS 

Consider  a  modification  of  the  NTC  collisional  sampling  approach  to  allow  for  collision  of  only  a  fraction  of  the 
number  of  real  particles  represented  by  each  computational  particle,  i.e.  a  fraction  of  the  statistical  ’weight’.  For 
every  collision  pair  selected,  the  collision  partner  weights  are  then  reduced  by  a  fraction  that  depends  on  the  collision 
probability  instead  of  performing  an  ’all  or  nothing’  collision  for  the  pair  depending  on  whether  a  random  number 
exceeds  a  probability  threshold,  the  defining  characteristic  of  ’Monte  Carlo’  methods.  The  fraction  of  real  particles 
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participating  in  the  collision  is  then  applied  to  newly  created  particles  such  that  the  result  of  the  collision  conserves 
the  total  weight. 

Although  more  conceptually  straightforward  than  the  sampling  only  a  subset  of  possible  collisions  and  randomly 
selecting  whether  they  occur  as  in  the  original  NTC  method,  this  approach  rapidly  increases  the  number  of  computa¬ 
tional  particles  in  the  simulation.  The  approach  therefore  requires  an  efficient  conservative  particle  merging  technique 
that  maintains  important  non-equilibrium  features  of  the  velocity  distribution.  In  this  work,  the  authors’  approach  for 
particle  merging  described  in  the  companion  paper,  Reference  [4],  is  used  to  counteract  this  growth  in  computational 
particles. 

With  this  approach,  the  number  of  collisions  to  perform  out  of  the  possible  N(N-l)/2  combinations  becomes  another 
tunable  parameter  for  the  simulation.  The  few  particles  that  occupy  the  tails  of  the  velocity  distribution  have  a 
disproportionately  high  impact  on  the  higher  order  moments  of  the  the  distribution  just  as  outliers  in  a  statistical 
distribution  impact  the  RMS  value  more  than  the  mean.  For  the  simulations  presented,  the  number  of  collisions  selected 
was  set  by  the  formula  shown  in  Equation  1 .  For  this 

N collision  =  ;jfA Jplog(Np)  (1) 

As  the  number  of  fractional  collisions  sampled  is  increased,  the  magnitude  of  the  change  of  moments  becomes 
smaller  for  any  given  timestep.  Rare  collisions  involving  the  distribution  tail  that  would  otherwise  only  occur  every 
few  timesteps  if  at  all  instead  scatter  a  small  fraction  of  the  probability  in  the  tails  every  timestep.  Even  though  the 
average  number  of  particles  scattered  out  of  the  tail  is  equal  across  many  timesteps,  it  is  hoped  that  the  fractional 
collisions  offer  a  possibility  for  smoother  and  more  tightly  coupled  collisional  evolution. 


BI-MAXWELLIAN  THERMALIZATION  RESULTS 

To  examine  the  performance  of  this  new  collision  sampling  technique,  consider  the  thermalization  of  a  bi-Maxwellian 
distribution  with  equal  y-  and  z-temperatures,  but  a  bimodal  x- temperature.  For  the  cases  presented,  the  thermal 
components  of  the  distributions  were  initialized  to  300K  in  y-  and  z-directions  and  150K  in  the  x-direction.  The 
centers  of  the  halves  of  the  bimodal  distribution  were  set  to  ±750  m/s  and  the  gas  used  was  Argon.  Collisions  were 
based  on  a  VSS  collision  cross  section  with  parameters  from  the  appendix  of  Reference  [1]. 

Figure  1  shows  that  the  rate  of  thermal  relaxation  for  10  separate  cases  using  the  identical  initial  particle  distributions 
for  the  original  NTC  method  and  the  new  fraction  collision  algorithm.  The  original  NTC  method  uses  800  particles 
for  this  test  case.  For  the  fractional  algorithm,  the  fraction  parameter  /  =  10-1/2  was  used  for  the  example  shown.  The 
merge  algorithm’s  target  number  of  particles  was  800,  but  the  actual  average  number  of  post  merge  macro  particles  was 
only  587.  This  number  was  the  result  of  a  balance  with  the  fraction  parameter  and  merge  aggressiveness  parameters. 
The  fractional  collisions  resulted  in  an  average  peak  of  2914  particles  prior  to  merging  back  down  to  587. 

Figure  2  shows  these  same  10  cases  where  the  mean  temperatures  and  RMS  of  the  deviation  from  the  mean  for 
each  timestep  were  calculated.  The  lines  and  error  bars  from  the  plots  are  these  means  and  ±RMS  values  respectively. 
Though  the  mean  temperatures  are  quite  consistent,  the  variance  of  the  temperature  components  in  time  match  those 
of  a  NTC  run  with  many  more  computational  particles. 

To  quantify  the  correlation  between  RMS  temperature  error  and  the  number  of  particles  at  various  stages  of  the 
fraction  collision  and  merge  process,  a  matrix  of  original  and  fractional  collision  runs  were  performed.  The  original 
collision  model  was  run  with  50,  100,  200,  400,  and  800  initial  particles.  The  fractional  collision  model  was  run 
using  the  original  800  particles  from  the  800  particle  original  case,  but  with  target  particle  numbers  of  50,  100,  200, 
400,  and  800.  This  array  of  target  particle  numbers  was  also  run  at  three  different  collision  fraction  parameters,  /, 
of  10_1,  10-1/2,  10°.  The  RMS  temperature  error  was  then  collected  for  each  of  these  cases  as  well  as  the  pre-  and 
post-collision  number  of  macro  particles. 

Figure  3  shows  the  effect  of  varying  the  number  of  particles  on  the  RMS  error.  The  variance  in  temperature  of 
the  new  method  appears  to  be  controlled  nearly  independently  from  the  pre-collision  number  of  particles.  Instead  it 
depends  on  the  peak  number  of  particles  resulting  from  fractional  collisions  prior  to  merging.  Because  the  number  of 
collision  partners  to  select  out  of  the  possible  N(N-l)/2  pairs  is  an  independent  parameter,  the  temperature  variance  can 
be  controlled  directly  at  the  expense  of  more  aggressive  particle  merging  at  the  end  of  the  collision  step  as  quantified 
by  the  /  parameter. 


FIGURE  1.  Thermal  relaxation  of  anisotropic  temperature  for  initially  bi-Maxwellian  temperature  distribution.  Plot  shows  10 
separate  cases  from  the  same  initial  particle  velocity  distribution.  Results  for  the  original  NTC  collisions  are  shown  on  the  left  and 
fractional  NTC  collisions  are  on  the  right. 


FIGURE  2.  Thermal  relaxation  of  anisotropic  temperature  for  initially  bi-Maxwellian  temperature  distribution.  Plot  shows  mean 
and  RMS  error-bars  of  10  separate  cases  from  the  same  initial  particle  velocity  distribution.  Results  for  the  original  NTC  collisions 
are  shown  on  the  left  and  fractional  NTC  collisions  are  on  the  right. 


COLLISIONAL  BEAMS  IN  POTENTIAL  WELL 

Having  shown  that  at  least  thermalization  rate  in  the  0D  sense  matches  between  the  original  and  fractional-NTC 
collision  algorithms,  it  is  important  to  determine  whether  these  the  fractional  collision  and  accompanying  merge  can 
reproduce  the  evolution  of  particles  that  are  also  moving  between  cells  spatially.  To  test  the  fractional  collision  in  ID 
rather  than  0D,  a  charged  particle  version  of  the  bi-Maxwellian  is  placed  within  a  1-D  electrostatic  potential  well  such 
that  the  two  oscillating  beams  interact  when  they  cross  at  the  bottom  of  the  well.  Though  the  particles  are  charged 
in  the  sense  that  an  external  force  accelerates  them  towards  the  spatial  center,  only  the  elastic  particle  collisions 
are  considered  and  modeled  with  a  VSS  collision  cross  section.  Consideration  of  effective-Coulomb  collisions  is 
considered  later. 

Figure  4  shows  that  the  variable  soft  sphere  (VSS)  thermalization  of  the  oscillating  beams  matches  well  between 
the  two  collision  sampling  techniques.  The  figure  also  shows  that  both  the  NTC  and  new  fractional-NTC  can  produce 
similar  number  density  evolution  despite  dramatically  different  numbers  of  computational  particles  per  cell. 

Unfortunately,  the  authors  know  of  no  exact  solution  for  this  crossed  beam  thermalization  problem.  Higher  fidelity 
simulations  using  traditionally  accepted  particle  methods  are  left  to  future  work  to  determine  whether  the  fractional 
collision  provides  any  advantage  over  traditional  methods. 
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FIGURE  3.  Temperature  RMS  error  relative  to  the  mean  temperature  on  average  for  all  timesteps  relative  to  the  number  of 
macro  particles.  The  left  plot  shows  the  pre-collision  average  number  of  macro-particles  as  the  x-coordinate  while  the  right  plot 
shows  the  same  error  levels  with  the  post-collision  average  number  of  macro-particles.  The  line  colors  depict  the  fractional-NTC 
aggressiveness  parameter. 


FRACTIONAL  MCC  IONIZATION  EVENTS 

For  a  second  test  of  fractional  collisions  with  merging,  a  electron  avalanche  in  a  spark  gap  gas  breakdown  is  considered 
as  it  was  in  the  companion  paper,  Reference  [4].  The  simulation  is  initialized  with  a  6kV  potential  drop  across  a 
5mm  gap.  The  evolution  of  electron  density  is  studied  in  a  partially  ionized  background  Argon  gas  which  is  assumed 
stationary  for  the  duration  of  the  simulation. 

The  inelastic  MCC  collision  model  previously  used  for  the  electron  avalanche  is  modified  to  remove  the  random 
nature  in  place  of  a  fractional  collision  as  it  was  with  the  NTC  elastic  collisions.  The  number  of  collisions  per  timestep 
is  calculated  exactly  as  it  was  in  Reference  [4].  Instead  of  selecting  collisions  that  exceed  a  probability  ratio  based  on 
that  collision  frequency,  the  number  of  colliding  particles  is  removed  from  the  original  particle  weight  and  two  new 
particles  are  created  following  the  formula  shown  in  Equation  2. 

Ar°  +  -F  ( 1  -  a)  (Ar0  +  e~rigina^j  +  a  (Ar+  +  e~attered  +  e~mitted)  (2) 

One  of  the  new  particles  represents  the  fraction  of  initial  particles  involved  in  the  collision.  The  velocity  of  this 
particle  is  reduced  by  the  ionization  energy  necessary  for  the  collision.  This  new  computational  particles  weight  is 
simply  the  number  of  collisions  per  timestep  used  to  reduce  the  original  particle’s  weight. 

The  second  of  the  new  particles  represent  the  electron  ejected  from  the  background  neutral  gas.  This  particle  again 
has  a  numerical  weight  equal  to  the  number  of  collisions  per  timestep,  but  its  velocity  is  zero. 

Though  this  collision  algorithm  triples  the  number  of  computational  particles  in  every  timestep,  the  merging 
algorithm  effectively  maintains  the  target  number  of  computational  particles  preventing  runaway  growth.  Figure  5 
compares  the  results  for  the  original  MCC,  merged  MCC,  and  merged  fractional-MCC  algorithms  for  the  spark  gap. 
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FIGURE  4.  Position  vs  Time  results  of  number  density  and  computational  particles  per  cell  in  thermalizing  oscillating  beams. 
The  left  figures  are  for  the  original  NTC  method  with  fixed  particle  weight.  The  right  figures  are  for  the  the  fractional-NTC  collisions 
along  with  octree  merge  algorithms. 


Two  key  differences  between  the  fractional  and  original  MCC  algorithms  are  apparent.  First,  the  fractional  MCC 
algorithm  has  less  ’streaks’  because  the  ionization  occurs  more  frequently  with  smaller  magnitude.  The  most  apparent 
streak  is  the  vertical  yellow  stripe  crossing  the  original  rarefaction  edge  at  approximately  4mm  and  0.2ns  in  the 
merged  MCC  example.  The  other  key  difference  is  that  the  fractional-MCC  results  become  noisier  in  general  above 
approximately  0.3ns.  This  is  likely  the  result  of  a  considerably  more  aggressive  merge.  With  the  fractional-MCC 
algorithm,  every  electron  above  the  threshold  energy  results  in  the  production  of  2  slower  particles  with  potentially  tiny 
weights  due  to  extremely  small  cross  sections  and  low  probability.  This  large  quantity  of  nearly  weightless  particles 
may  be  overwhelming  the  octree-merge.  A  higher  fidelity  reference  case  is  needed  to  further  quantify  ’correctness’  in 
the  absence  of  an  exact  analytical  solution  for  the  ionizing  spark. 

Though  the  electron  density  evolution  looks  similar  between  the  MCC  and  fractional-MCC  collision  methods,  one 
area  in  which  the  fractional-MCC  method  particularly  excelled  was  in  the  smoothness  of  the  positive  ion  profile. 
Because  the  fractional  MCC  collision  event  creates  particles  with  only  a  fraction  of  the  weight  of  the  original,  the 
ion  density  is  incremented  in  a  smooth  fashion  between  cells  and  timesteps.  However,  again  because  no  analytical  or 
very  high  fidelity  solution  exists,  it  is  impossible  to  determine  if  these  results  are  ’better’  or  not  with  the  current  tools. 
Again,  higher  fidelity  baseline  simulations  are  needed  to  determine  whether  the  fractional  collision  model  is  actually 
more  predictive. 


FUTURE  EXTENSION  AND  CONCLUDING  REMARKS 

In  addition  to  studying  fractional  elastic  and  inelastic  ionizing  collisions,  a  model  for  an  effective  Coulomb  collision 
is  also  currently  under  consideration.  Without  such  a  model,  the  populations  of  free  streaming  and  post-collision  or 
post-ejection  electrons  have  no  mechanism  to  thermalize. 

Unfortunately,  the  cross  section  for  Coulomb  collisions  diverges  without  a  cutoff  to  the  integral  at  the  Debye  radius. 
This  cutoff  is  ill-defined  in  a  highly  non-equilibrium  sense.  Short  of  modeling  the  full  N-body  interaction,  or  at  least 
the  multi-pole  expansion  thereof,  a  simplified  collision  model  is  needed. 
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FIGURE  5.  Position  vs  Time  diagram  of  electron  density  and  computational  particles  per  cell  in  an  ionizing  spark  gap  using 
MCC  collision  model.  Electrons  are  pushed  from  the  cathode  on  the  left  to  anode  on  the  right  of  the  figures  due  to  the  6kV  applied 
potential.  The  left  figure  shows  the  original  MCC  collision  model  with  fixed  particle  weight  and  the  middle  figure  shows  the  original 
MCC  collision  model  for  the  octree  merge  algorithm.  The  right  figure  shows  the  fractional-MCC  model  with  octree  merge. 


Such  a  simplified  approximate  model  based  on  ZeTdovich’s  approach  to  estimating  electron  conductivity  [5]  is 
considered.  Equation  3  shows  this  effective  cross  section  for  Coulomb  collisions  from  Zel’dovich,  though  the  leading 
constant  is  given  without  justification.  The  model  estimates  an  effective  cross  section  and  collision  frequency  for  which 
the  accumulation  of  many  small  angle  Coulomb  collisions  results  in  a  momentum  change  of  the  order  of  the  original 
velocity.  The  Coulomb  logarithm,  InA ,  accounts  for  the  cutoff  to  bound  the  collision  integral  for  ’far’  interactions. 
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The  same  basic  approach  was  used  by  Kapper  and  Cambier  [6]  to  estimate  the  electron-ion  thermalization  in 
their  two-fluid  Argon  shock  ionization  collisional  radiative  studies.  As  these  results  compared  well  to  Schlieren 
and  interferogram  experimental  results,  applying  this  effective  approximate  Coulomb  collision  appears  to  be  a  good 
starting  point  for  electron  thermalization  and  conductivity  in  PIC  codes.  The  expectation  would  be  that  the  electron 
conductivity  would  the  at  least  be  consistent  between  particle  and  fluid  models  near  equilibrium  and  would  provide 
initial  insight  into  the  potential  breakdown  when  moving  away  from  equilibrium. 

For  an  initial  study  of  self-thermalization  of  electrons,  the  cross  section  of  Equation  3  was  modified  by  replacing  the 
temperature  with  relative  velocity  as  kT  =  \ mv 2.  Though  only  approximately  correct,  the  goal  was  to  at  least  achieve 
the  correct  near-equilibrium  thermalization  rate  when  the  temperatures  are  well  defined. 

Several  OD  electron  thermalization  cases  were  performed,  but  relaxation  rates  appeared  to  be  highly  dependant  on 
time  step  and  merge  parameters.  Though  further  investigation  is  needed, 

it  is  thought  that  this  result  may  be  due  to  the  oc  T-2  dependence  of  the  cross  section  on  temperature  in  conjunction 
to  the  particle  merge.  The  hard  sphere  model  has  a  temperature  dependence  related  to  T1/2  and  the  VHS  and  VSS 
models  vary  between  T1/2  and  T1.  The  poor  performance  of  the  pseudo-Coulomb  collision  might  be  a  result  of  the 
merge  only  conserving  moments  through  v2  while  the  Coulomb  cross  section  depends  on  v4.  The  high  order  velocity 
dependence  means  that  the  effective  Coulomb  relaxation  is  highly  dependent  on  the  shape  of  the  tail  of  the  distribution 


that  is  not  necessarily  preserved  by  the  merge.  Hopefully  this  effect  can  be  mitigated  in  the  future  with  a  merge  that 
conserves  higher  moments  as  initially  outlined  in  the  companion  merge  paper[4] . 
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o  Ionizing  Breakdown  in  6kV 
Potential 

o  Electron  Flow:  Cathode  — >>  Anode 
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Fractional-MCC  Gas  Breakdown 
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o  Ionizing  Breakdown  in  6kV 
Potential 

o  Electron  Flow:  Cathode  — >>  Anode 
o  Inelastic  Fractional-MCC  Version 
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o  Ionizing  Breakdown  in  6kV 
Potential 

o  Electron  Flow:  Cathode  — >>  Anode 
©  Inelastic  Fractional-MCC  Version 
©  Fess  Streaks  Early 
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o  Ionizing  Breakdown  in  6kV 
Potential 

o  Electron  Flow:  Cathode  — >>  Anode 
©  Inelastic  Fractional-MCC  Version 
©  Fess  Streaks  Early 
©  More  Noise  Eater 
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Pseudo-Coulomb  Collisions 
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Difficulty  with  Coulomb  Collisions: 

o  Cross  Section  Diverges 

o  Cannot  be  Fully  Represented  Without 
N-body  Interactions 
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N-body  Interactions 

o  Surrogate  needed  for  Coulomb 
Thermalization 
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Cross  Section  Diverges 

Cannot  be  Fully  Represented  Without 
N-body  Interactions 

Surrogate  needed  for  Coulomb 
Thermalization 

Zel’dovich  describes  Combined 
Collisions  «  to  Original  Momentum 

Non-Equilibrium  Cross  Section  to 
Reproduce  Rate  Attempted 

Early  Results  AT/Merge  Dependent 
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Pseudo-Coulomb  Collisions 


Difficulty  with  Coulomb  Collisions: 


o  Cross  Section  Diverges 

Ai 

o  Cannot  be  Fully  Represented  Without 
N-body  Interactions 

o  Surrogate  needed  for  Coulomb 
Thermalization 

o  Zel’dovich  describes  Combined 

Collisions  «  to  Original  Momentum 

N 

o  Non-Equilibrium  Cross  Section  to 
Reproduce  Rate  Attempted 

o  Early  Results  AT/Merge  Dependent 

o  Cross-Section  oc  T~ 2  «  v-4 

Must  Conserve  Higher  Order  Moments? 

/T  ^  HT  ^  " 

°  ~  9^ 
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Collision-Radiative  Shock  Ionization 
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Extending  Applicability  from  2-Fluid: 

o  Cellular  Structures  and  Electron  Avalanche 
Phenomena  found  in  Shock-Ionization  of 
Argon 

o  2-Fluid  Model  Ignores 

Rarefaction/Non-Equilibrium  Behavior 

o  Particle  Code  would  Enable  Rarefied  Studies 

©  Future  Hybridization  of  2-Fluid/Particle 
would  Enable  Engine  — >  Plume  modeling  of 
Spacecraft  Electric  Propulsion 


Kapper  &  Cambier,  J.  Appl.  Phys.  109, 
(2011) 
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